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ABSTRACT 

Compact plasmas, that exist near black-hole candidates and in gamma ray bursts 
sources, commonly exhibit self-organized non-linear behavior. A model that simulates 
the non-linear behavior of a compact radiative plasma is constructed directly from the 
observed luminosity and variability. The simulation shows that such a plasma self- 
organizes, and that the degree of non-linearity as well as the slope of the power density 
spectrum increases with compactness. The simulation is based on a cellular automaton 
table that includes the properties of the hot (relativistic) plasma, and the magnitude 
of the energy perturbations. The plasma cools, or heats up, depending on whether it 
releases more or less energy than that of a single perturbation. The energy released 
depends on the plasma density and temperature, and the energy of the perturbations. 
Strong perturbations may cool the previously heated plasma through shocks and/or 
pair creation. 

New observations of some active galactic nuclei and gamma ray bursts are consistent 
with the simulation. 
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1. INTRODUCTION 

Large amplitude X-ray and 7-ray variabil- 
ity from active galactic nuclei (AGNs), black 
hole binaries (BHBs), and gamma ray bursts 
(GRBs) are thought to be the signatures of hot 
(relativistic) plasma (Blandford 1990, hereafter 
B90, and references therein, Mushotzky, Done 
& Pounds 1993, hereafter MDP93, Fishman & 
Meegan 1995, hereafter FM95). Many of these 
sources exhibit rapid large amplitude variabil- 
ity (MDP93, Green, McHardy & Lehto 1993, 
hereafter GML93). The combination of rapid 
large amplitude variability and high emissivity 
is thought to be evidence for compact sources 
(Fabian 1992, B90, MDP93, FM95). 

Roughly half of AGNs and BHBs exhibit l/f a 
decline in their power density spectrum (PDS), 
where 1 < a < 2 (Press 1978, MDP93, McHardy 

1988, GML93, Ulrich et al. 1997). Some GRBs 
also exhibit non-linear oscillations superposed on 
the general power law decay (Meredith, Ryan & 
Young 1995, hereafter MRY95). Recently sev- 
eral time series analysis methods were devel- 
oped to distinguish between PDS with linear 
origins and PDS with non-linear origins (e.g., 
Kaplan & Glass 1995, hereafter KG95, Vio et 
al. 1992, hereafter V92, Scargle 1990). By 
using some of these methods evidence of non- 
linearity in some AGNs and GRBs was found 
(V92, Boiler et. al. 1997, Leighly & Obrian 1997, 
MRY95, Yuan et al. 1996). Some of these light 
curves were also found to be self-similar (Scar- 
gle, Steiman-Cameron & Young 1993, hereafter 
SSY93, McHardy & Czerny 1993). 

Several authors have borrowed cellular au- 
tomaton (CA) models from other scientific fields, 
to simulate the non-linear behavior of AGNs, 
BHBs and GRBs. CAs are tables of rules that 
describe how to model a non-linear system with 
identical elements. CAs are commonly used 
when the partial differential equations for a non- 
linear system are not easily solvable (Jackson 

1989, hereafter J89, KG95). Using a CA Mi- 



neshige, Takeuchi & Nishimori (1994, hereafter 
MTN94), assumed that the accreting material 
in BHBs is in the form of an accretion disk, 
and that the avalanches are analogues to those 
in the self organized criticality (SOC) CA for 
sand piles (Bak, Tang &; Wicscnfcld 1988, here- 
after BTW88). SSY93 assumed that the mate- 
rial in BHBs is in the form of a ring, as in the 
dripping hand model of Crutchfield &; Kaneko 
(1988). Stern k Svensson (1997, hereafter SS97) 
suggested that a pulse-avalanche CA may be ap- 
propriate for an expanding "fireball" in GRBs if 
there is magnetic turbulence due to an unknown 
instability. These models are motivated by simi- 
larity to simulations from other fields, and by the 
assumed geometry of some astrophysical point 
sources. The aim of this Letter is to construct 
an independent CA directly from the accepted 
physical conditions in the plasma. 

The physical conditions in compact plasma 
with perturbations and the CA table of rules are 
described in section §2. In section §3 it is shown 
that the plasma indeed evolves to a self-organized 
critical state for a wide range of emission rates 
and that the emerging power-law spectrum is of 
the non-linear l/f a variety. In section §4, the 
last section, the expected observational ramifica- 
tions and conclusions are discussed. 

2. CA MODEL FOR COMPACT PLAS- 
MAS WITH PERTURBATIONS 

In many AGNs, BHBs and GRBs there are 
compact plasmas with nonlinear processes. The 
energy lost to emission or added to the plasma 
on a doubling time scale is a significant fraction 
of the total energy (Sivron 1995, hereafter S95, 
Piran 1995, hereafter P95). The introduction 
of a large energy density perturbation into the 
plasma should result in an increase of emission 
if the energy is efficiently radiated. For example, 
emission will rise due to a pair-runaway process 
in which the increased number of leptons results 
in efficient Compton cooling. 



2 



Post-shock pair cascades are readily produced 
in compact sources with strong pertrubations 
(Sivron, Caditz & Tsuruta, hereafter SCT96). 
Strong perturbations, in which the perturbation 
excess density is of the order of the average den- 
sity in the plasma are known to form shocks 
(Landau & Lifshitz 1987). This holds true for rel- 
ativistic fluids and hot collisional plasma (Taub 
1949, Iwamoto 1989, hereafter 189). In many 
situations the time scale for perturbation steep- 
ening to shocks is small. This is known to be 
true for compact accreting sources (Papaloizou & 
Pringle 1984, Narayan 1991), and is conjectured 
for GRBs (P95). In such cases one may model all 
strong super-sonic perturbations as effective ra- 
diative shocks that can dominate the light curve 
(SCT96). The same strong perturbations may 
result in very little radiation, when the plasma 
temperature is too low for pairs to be effectively 
created. In such a case the perturbation energy 
would heat up the plasma. All the perturbations 
are nonlinear in the sense that a strong pertur- 
bation significantly changes the plasma temper- 
ature and speed of sound, which then determine 
the output due to the next perturbation. 

In the simulations in this paper it is assumed 
that non-linear strong perturbations are pro- 
duced in the source. The input for the simula- 
tions is a perturbation moving with input veloc- 
ity Ui , and the output is the luminosity Lj , where 
the number of intervals is N, and the running in- 
dex is 1, 2, i, N. The simulations follow the 
CA rules of table 1. The simulation is non-linear 
in the sense that parameters in cells i + i + 2 
etc., depend on the output Lj. The CA rules and 
the parameters in it are described in the next 
three paragraphs. 

For demonstrative purposes numbers which 
are appropriate for a typical source, the active 
nucleus of a Seyfert I galaxy, are used. The 
emitting source of size X = 5Rs c h = 1-5 x 
10 13 cm, appropriate for a central black hole of 
mass M = 10 7 Mq, and the accretion rates rel- 
ative to the Eddington accretion rate are m = 



0.1,0.5,1.0,1.5. Here R Sc h = 2GM/c 2 is the 
Schwarzchild radius of the putative central black 
hole with mass M, rh = M/Mem-i M is the ac- 
cretion rate, M.Edd = ^GMmn / (car) is the Ed- 
dington accretion rate, and run and ot are the 
hydrogen mass and the Thompson cross section 
respectively (B90). 

The parameters for the simulation are the fol- 
lowing: The time interval is A, and A = 1 second 
for the case of Seyfert Is. The total mass of the 
plasma ra = 3 x 10 24 gram, was chosen so that, 
using the plasma deflection length with electron 
temperature ~ 3 x 10 9 K, the plasma deflection 
length is smaller than X, and the plasma is colli- 
sional (189). The total energy of the plasma, not 
including rest mass, is Ei and the total temper- 
ature of the plasma is Tj = Ei/(1.5riik), where 
Ui = 3m/(m#47rX 3 ), k is the Bolzman constant 
and ran is the average mass of an atom. The 
mass of each perturbation is rrii = raA, and 
the kinetic energy of a perturbation is (SE)i = 
[0, (5E) max ] with equal probability^. (5E) max is 
the Keplerian energy of perturbation with mass 
rrii at 5Rs c h- Using the usual special relativis- 
tic expression, the input speed of each perturba- 
tion is Ui = Cyjl — (rriiC 2 /(5E)i) 2 . The output 
luminosity is Li , and the efficiency of converting 
gravitational to radiative energy is ef]. The initial 
conditions were set at multiples of E\ = 0.05mc 2 
(see results, section 3). Energy perturbations 
with Gaussian distributions were also tried. The 
results were generally the same (see section 3). 

The effective shocks are selected in the fol- 
lowing way: The parameters Ti,rii,Ui are sent 
to a subroutine that uses the Ranking-Hugoniot 
relations for a hot collisional plasma to find the 
post-shock conditions (189, SCT96). If these con- 
ditions are sufficient for pair-cascades, and if the 
perturbation moves faster than (c s )j, the shock 

2 The distribution is the result of virial radiative pressure 
and gravity that produce super and sub-Keplerian speeds 
for perturbation of large enough size (S95, SCT96). 

3 In the case of accretion sources e = 0.06 for a 
Schwarzchild black hole. 
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is considered effective (Svensson 1982, Svensson 
1984, BS91, SCT96). The speed of sound in 
the plasma, (c s )i depends on Tj, and is typically 
slightly larger than c/3. There are two alternat- 
ing CAs: If there are no effective shocks the CA 
is the second row in the table. If there is an effec- 
tive shock the third row is the CA. The number 
of columns in that row is J = (X/Ui)/A + 1, 
an integer that determines the number of subse- 
quent cells still affected by effective shocks. In 
table 1 the width of row B was selected to be 
J = 3 for demonstrative purposes. 

The following scenario demonstrates how the 
non-linearity of the simulation: If initially there 
is no effective shock the energy Ei+\, is prob- 
ably smaller than E{. In such cases LjA/e is 
much smaller than (5E)i+\, because the small 
non-shock emission is diffusion - dominated (see 
table 1). As a result the energy and tempera- 
ture of the plasma in the subsequent intervals 
increase, as does the speed of sound. In the fol- 
lowing interval the probability of exceeding the 
speed of sound is therefore smaller. If, on the 
other hand, there is an effective shock, the en- 
ergy decreases according to columns 2 — 4 in the 
third row, lowering the speed of sound, and in- 
creasing the probability of shocking the plasma 
in the next interval. 

3. RESULTS 

The light curve, PDS, phase, and auto - cor- 
relation function (ACF) of the output from the 
model in Figs la, lb and lc resemble the re- 
sults of BTW88, MTN94, SSY93 and SS97. As 
in BTW88 self organization is achieved indepen- 
dently of initial conditions, because the system 
quickly evolves to a state in which perturbations 
of above average energy shock the plasma. Self 
organization is not obtained, however, for accre- 
tion rates lower than m = 0.1 or larger than 
m = 1.0. At m < 0.1 the plasma temperature 
grows without bound because there are too few 



large perturbations at supersonic speeds^. For 
m > 1.0 our model was not stable. 

In Fig la the light curve corresponding with 
the higher m takes longer to reach the organized 
criticality state because while increasing its tem- 
perature the denser plasma loses more energy 
through small shocks. The minimum doubling 
time scale is smaller for the higher ra objects 
because we increased rh by lowering the size of 
the plasma. The shock therefore extracts energy 
from the plasma over a shorter period, resulting 
in stronger flares. For accreting sources this cor- 
responds with a smaller compact object. 

The PDS exponent in Fig lb is a = 0.85 ±0.03 
for m = 1 and a = 0.79 ± 0.03 for m = 0.5 
(fit not shown). In both cases the PDS in- 
cludes emission from the time after the criti- 
cal state has been established, and the subse- 
quent 900 seconds. Energy perturbations with 
Gaussian distribution of width 0A5E max yielded 
a = 0.84 ± 0.03 and 0.79 ± 0.03. As expected, 
for increasing m small high frequency shocks are 
increasingly suppressed as the temperature and 
associated speed of sound increase just before the 
more frequent large shocks. With decreasing m 
we get a ~ 0, because as the time interval be- 
tween perturbations grows the system responds 
linearly to the random perturbations. These re- 
sults are similar to those in BTW88. In BTW88 
a flat PDS is changed to l/f a due to suppression 
of short scale avalanches in domains of increas- 
ing sizes that were flattened by larger avalanches. 
In the model presented here the l/f a is due to 
the suppression of weak radiative shocks just be- 
fore stronger shocks. Here the "critical slope" of 
BTW88, the speed of sound, depends on the pre- 
vious temperature in the plasma. This model is 
therefore analogous to a one dimensional BTW88 
sand pile with critical angle that depends on the 
speed and location of previous avalanches. 

4 In such cases modest magnetic fields are needed for cool- 
ing with super Alfvenic perturbations and reconnection 
events. 
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Correlation on short time scales (in the first 
few seconds and at around 50-80 seconds) for 
rh = 1 can be seen in the ACF in Fig lc. The 
enhanced ACF at low time scales makes sense be- 
cause of the anticipated correlation of perturba- 
tions. The peaks are due to the total sum of dif- 
ferent average delays due to different conditions 
in the plasma. The effect of the initial pertur- 
bation is lost over time scales t > X/{Aii{c s )i), 
where A4j = Ui/(c s )i is the Mach number. For 
rh = 0.5 there is less overall correlation, but more 
correlation at times ~ 80 seconds (not shown). 

Another method by which the non-linear de- 
pendence is demonstrated is shown in Fig 2a in 
which the phase space diagram of the outputs L jy 
is compared with L/v+i- With no correlation the 
path should randomly fill the correlation space, 
as seen in Fig 2b for a 100 second delay between 
Ljy with Lat + ioo- Despite of the apparent cor- 
relation for a 70 second delay the phase space 
diagram does not show obvious structure (not 
shown) . 

4. DISCUSSION AND CONCLUDING 
REMARKS 

The model presented here shows that l/f a 
non-linear PDS can be created without relying 
on a specific geometry. One problem with this 
model is that, in accreting sources, the pre- 
dicted a is too low. This result can be cor- 
rected by adding parametric dimensions. Adding 
parametric dimensions usually results in larger a 
(BTW88). We are currently working on a simu- 
lation with added parameters, one of which is 
related to the 2D angular momentum transfer 
- an essential component for accreting sources 
(Sivron & Leighly 1998, hereafter SL98, Sivron & 
Tsuruta 1993, hereafter ST93). For low rh and 
flat geometry the added parameter is expected to 
yield results similar to those of MTN94, because 
the random input and output for each disk cell 
are correlated. A parameter related to the pro- 
files of radiation events from each shock is ex- 
pected to yield a smaller a, due to relativistic 



effects (SL98). This parameter is analogous to 
the time profile of shots in Takeuchi, Mineshige 
& Negoro 1995. Another problem with the model 
is that, contrary to observations, there is no cor- 
relation on a time-scale longer than 10 -3 seconds 
for BHBs (Negoro et al. 1994). This problem 
will be corrected with more dimensions because 
rh will decrease with the less effective angular 
momentum transfer associated with small non- 
shocking perturbations, making the subsequent 
shocks less effective (SL98). 

For accreting sources the simulation yields av- 
erage output temperatures that are lower for high 
rh. This is because of an increase in emission ef- 
ficiency with increasing rh. However, when rh 
is high the cooling of the post-shock plasma is 
effective, and a large portion of the post-shock 
material is cooled to a "cold phase" of temper- 
ature ~ 10 6 K. The effect of such matter on 
observations is an enhancement of the soft X-ray 
emission. (Guilbert & Rees 1988, Celotti Fabian 
& Rees 1991, Sivron & Tsuruta 1993, SCT96, 
Kuncic, Celotti & Rees 1997). This result is al- 
ready consistent with observations of narrow line 
Seyfert galaxies type I (NLSIs) and "regular" 
Seyfert Is. NLSIs that have the same emission 
rate of regular Seyferts have steeper X-ray spec- 
trum, a slightly larger a, and probably smaller 
black hole and larger rh (Leighly 1997, SL98). 

The correlation of rh, a and enhanced ther- 
mal emission is general. The model therefore 
predicts that energy density perturbations in the 
initial fireball of GRBs can produce non-linear 
temporal variations. These variations can then 
give rise to winds with varying Lorenz factor T 
(S95, SL98). The emission can be the result of 
low T winds loaded with baryonic matter over- 
taking the high T initial winds that slow down 
as the fireball sweeps up external matter (Rees 
& Meszaros 1997). The non-linear variations in 
this scenario are frozen in, and reveal themselves 
when the initial fireball expands and becomes 
optically thin. This scenario is consistent with 
GRBs observed light curves. 
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Fig. 1. — a. The upper light curve corresponds 
with the higher m = 1, and takes longer to reach 
the organized criticality state, but the minimum 
doubling time scale is smaller. The lower curve 
corresponds with fa = 0.5. b. In the lower box 
the solid curve represents the PDS. The dashed 
curve is the fit to the PDS in the case m = 1, 
with a = 0.86 ± 0.03. In the upper box the 
solid curve is the phase of the PDS. c. The auto- 
correlation function for rh = 1 (bold) is com- 
pared with the auto-correlation function for the 
more jagged random input, which has a typical 
exponential decay. There is an enhancement in 
correlation in the first few seconds and at around 
50-80 seconds. 



Fig. 2. — a. Phase space diagram with 1 second 
delay shows a very clear sign of non-linear depen- 
dence. The numbers on all axes are in arbitrary 
units. The order of points is such that for a point 
N on the lower accumulation line the point iV+ 1 
is on the upper accumulation line. b. The same 
phase relation with 100 second delay shows no 
dependence at that time interval, corresponding 
with loss of information over times larger than 
X/U N . 



This 2-column preprint was prepared with the AAS IATgX 
macros v4.0. 
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Table l:The CA Rules 



Conditions for shock 


L at cell i 


Effects of Ei 


Effects of Ei 


on 


(T i = T(E i ),cs i = cs(T i )) 




on cell i + 1 


on cell % + 2 


i + 3 




t = iA 


t = (i + l)A 


t = (z + 2)A 


t = (i + 3)A 


IF Ui < cs** 


U = dip 


E l+1 =Ei + (SE)i 


no 


no 


no effective shock 




-Li * A/e 


effect 


effect 


If Ui > cs** 


Li = 0.5eEi 


Ei + {8E)i/J- 




E i+2 + (^)/J- 


shock is effective 


/(«/ + l) ft 


Li * A/e 


Lj * A/e 


Lj * A/e 



* The approximate expressions are: Tj ~ 2Ei/3k, 
csi ~ y/2kTi/m. (The fully relativistic expres- 
sions are from 189.) 



** More accurately, effective shocks are only 
those shocks in which Ui > csi and copious pair 
production is achieved (see SCT 96). 
tThe results of White & Lightman (1989) for 
Comptonized bremsstrahlung emission in hot two 
- temperature plasmas for a given rh are used for 
the diffusive emission dif. Hot two-temperature 
plasmas are expected for rh ~ 1 (Rees et. al. 
1982, Narayan & Yi 1994). 

ft The 0.5 parameter: due to assumption that half 
of the energy in the plasma at time ti is radiated 
away following a shock (see text). 
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